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Abstract. A stochastic theory for the toppling activity in sandpile models is 
developed, based on a simple mean-field assumption about the toppling process. The 
theory describes the process as an anti-persistent Gaussian walk, where the diffusion 
coefficient is proportional to the activity. It is formulated as a generalization of the Ito 
stochastic differential equation with an anti-persistent fractional Gaussian noise source. 
An essential element of the theory is re-scaling to obtain a proper thermodynamic limit, 
and it captures all temporal features of the toppling process obtained by numerical 
simulation of the Bak-Tang-Wiesenfeld sandpile in this limit. 
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1. Introduction 

The existence of self-organized critical dynamics in complex systems has traditionally 
been demonstrated through numerical simulation of certain classes of cellular automata 
referred to as sandpile models [lj. Non-linear, spatio-temporal dynamics is always 
essential for the emergence of SOC behavior, but the details of this dynamics for a 
specific natural system is often poorly understood and/or not accessible to observation. 
In many cases the information available is in the form of time-series of spatially averaged 
data like stock-price indices, geomagnetic indices, or global temperature data. For 
scientists who deal with such data a natural question to ask is: are there specific 
signatures of SOC dynamics that can be detected from such data? 

In this letter we shall report some results which provide a partial answer to such 
a question. Some important statistical features of the toppling activity are common to 
most weakly driven sandpile models described in the literature, and these are used to 
formulate a stochastic model for the toppling activity signal. A benchmark case against 
which our results are tested, is a numerical study of the Bak-Tang-Wiesenfeld (BTW) 
sandpile [2j . A crucial step in our work is a re-scaling of the dynamical variables which 
allows a natural passage to the thermodynamic (continuum) limit. We demonstrate 
that this leads to new results concerning SOC scaling laws. We find that the probability 
density function (pdf) for the toppling activity is a stretched exponential or close to 
the Bramwell-Holdsworth-Pinton distribution [3j, depending on whether the sandpile 
is so slowly driven that avalanches are well separated, or it is driven so hard that 
several avalanches run simultaneously. The pdf for avalanche durations is unique in 
the thermodynamic limit, but is not a power law, unless we redefine the meaning of an 
avalanche to be the activity burst between successive times for which the activity rises 
above a positive threshold. Implementing such a threshold yields an exponent for the 
avalanche duration pdf of 1.63, in agreement with pE], but in contradiction to [7J. It 
also gives power-law quiet-time statistics as in [1] and thus refuting the claim in [8] that 
SOC implies power-law distributed avalanche durations, but Poisson-distributed quiet 
times. 

The sandpile models considered in this short paper deal with ad> 2-dimensional 
lattice of N d sites each of which are occupied by a certain integer number of quanta 
which we conveniently can think of as sand grains. The dynamics on the lattice is 
given by a toppling rule which implies that if the number of grains on a site exceeds a 
prescribed threshold, the grains on that site are distributed to its nearest neighbors. If 
the occupation number of some of these neighbors exceed the toppling threshold these 
sites will topple in the next time step, and the dynamics continues as an avalanche until 
all sites are stable. The details of this toppling rule can vary, but a useful theory for a 
broad class of natural phenomena should not be very sensitive to such detail. 

In natural systems the SOC dynamics is usually driven by some weak random 
external forcing. In sandpile models this can be modeled by dropping of sand grains at 
randomly selected sites at widely separated times. In numerical algorithms this is often 
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Figure 1. a): A realization of the toppling activity xn(£) in the BTW sandpile. b): 
The increments Ace at (t) = X]y(t + 1) — Xw(t) of the trace in (a), showing that Ax^lt) 
is large when Xjv(i) is large, c): Conditional pdfs of xn + Axn for xn — 10, 20, 30 
respectively, d) The conditional mean and variance of Axn versus xn- 

done by dropping sand grains only at those times when no avalanche is running. This 
ensures that the drive does not interfere with the avalanching process. Usually it will 
then only take a few time steps from one avalanche has stopped until a new starts, so for 
a large system the quiet times between avalanches will appear insignificant compared 
to their durations. 

A more physical drive would be to drop sand also during avalanches. If the dropping 
rate is slower than the typical duration of a system-size avalanche the drive would still 
not interfere with the avalanche dynamics, but the quiet times would depend on the 
statistical distribution of dropping times, which is typically a Poisson distribution. In 
many natural systems, however, avalanching occurs all the time, corresponding to a 
higher driving rate. In such cases, and also because there will always be noise in time- 
series data, we cannot identify the start and termination of an avalanche from a zero 
condition of our observable. In practice we have to define avalanches as bursts in the 
time series identified by a threshold on the signal [I]. In a sandpile simulation such 
bursts are correlated and therefore the quiet times between the bursts are power-law 
distributed even if the dropping of sand grains is chosen to be a Poisson process. Hence 
if focus is on modeling features that can be detected in observational data we shall think 
of avalanches as activity bursts starting and terminating at a non-zero threshold value. 
Moreover, one of the main results of this work is that power-law shape of the pdf for 
avalanche duration is true only if one defines avalanches in this way. 

2. The stochastic model 

We shall assume that the lattice has linear extent L = 1 with N d sites, so the 
thermodynamic limit N — > oo can be thought of as a continuum limit. The sandpile 
evolves in discrete time steps labeled by k — 1, 2, 3 . . ., and the number of sites whose 
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occupation number exceeds the toppling threshold at time k is called the toppling 
activity x^(k). The toppling increment is Axx(k) = f x^{k + 1) — xjf(k). Let us define 
two active sites as dynamically connected if they have at least one common nearest 
neighbor, and define a connected cluster as a collection of active sites which are linked 
trough such connections. From numerical simulations of sandpiles we observe that such 
clusters never consist of more than a few elements and that the instantaneous number 
of clusters increases in proportion to xjy. This implies that at each time k we can 
label the clusters by i = 1, . . . , cx^ik), where c < 1 is a constant depending on the 
specific toppling rule and the dimension d of the sandpile. We can then decompose 
the increment Axw{k) into a sum of local increment contributions ^N,i(k) produced by 
each of the clusters, i.e. Ax^(k) = Y^Li ^N,i(k). We think of the local increment 
contributions as random variables which take values in a finite sample space. Indeed, if 
each cluster i only consists of a single overcritical site, then takes values in the set 
{-l,0,...,2d-l}. 

As a first step to a stochastic model we make a mean- field assumption [TOl [TT] . 
which impiles that £jv,i(fc) an d £,N,j(k) are statistically independent for i ^ j. Then the 
central limit theorem states that in the limit N — > oo, XN{k) — > oo the conditional 
probability density P[AxN(k)\x]sr(k)] of an increment Axx(k), given x^(k), is Gaussian 
with variance a 2 X]y(k), where a 2 = c 2 ''(E[^% { \x n] — (E[^ N<i \xiy}) 2 ). This has been 
verified numerically in the two-dimensional BTW-model as shown in Fig. [TJ The figure 
demonstrates the need to introduce a conditional probability: The conditional variance 
of the increments is proportional to xn and the conditional mean is not zero. 

In fact, numerical simulations show that the the conditional mean increment, 
E[Axn\xn], is positive for small xn, reflecting the natural tendency for the activity 
to grow when it is small. On the other hand the mean increment decays exponentially 
to zero for moderate x^, and becomes negative when x^ is comparable to the activity of 
a system-size avalanche, reflecting the limiting influence of the finite system size. These 
effects will be incorporated as a drift-term correction to the model, but for now we 
consider for simplicity of argument a Gaussian process with non-stationary increments 
and no drift term: 



Ax N (k) = a \/x N {k) w{k) , (1) 

where w(k) is a stationary Gaussian stochastic process with unit variance. From the 
numerical sandpile data (see Fig. [TJ we observe that the normalized toppling process 

def „ /N Ax N (k') 



k'=0 k'=0 °V X N(k') 

has the characteristics of a fractional Brownian walk with Hurst exponent H w 0.37 
on time scales shorter than the characteristic growth time for a system-size avalanche, 
consistent with a power spectrum which scales like /~ 174 . Thus we model the normalized 
increment process as w(k) = Wjj{k + 1) — Wn{k), where Wu(k) is a fractional Brownian 
walk with Hurst exponent H. For the transition to the thermodynamic limit, where 
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time will become a continuous variable, we can think about Wjj{k) as the result of a 
discrete sampling of the (continuous-time) fractional Brownian motion (ffim) Wn{t)- 
This process has the property (| Wn(t + r) — Wnit) | 2 ) = r 2H . We now have a stochastic 
difference equation 

Ax N (k) = a y/x N (k) (W H (k + 1) - W H (k)). (2) 

Numerical simulations show that x^ ~ N Dl where < D\ < d can be interpreted 
as a fractal dimension of the set of active sites imbedded in the ^-dimensional lattice 
space. This property is used to re-scale xjy(k) such that it has a well-defined limit 
as iV — > oo. We also have to re-scale the time variable by letting t = kAt, where 
At = N~ D2 . The value of D 2 will become apparent if we define the normalized activity 
variable X^(t) = N~ Dl XN{t/At), such that the corresponding increment becomes 

AX N (t) = N HD ^ a ^X^)AW H {t) , (3) 

where AWn{t) = Wjjit + At) — Wuit)- A well-defined thermodynamic limit N — > oo 
requires D 2 = Di/2H, for which Eq. ([3]), by introduction of the limit function 
X(t) = limjv^oo Xjv(t), reduces to the stochastic differential equation 

dX(t) = f(X) dt + a y/X(t) dW H (t), (4) 

where we have heuristically added a drift term f(X) dt to account for the non-zero 
mean of the conditional increment. We take f(X) to be an exponentially decaying 
function based on the numerical results from the sandpile. In the 2-dimensional BTW 
model we find that D\ pa 0.86 and hence D 2 = 1.16. This defines re-scaled coordinates 
X N = x N /N om and t N = k/N 1AG . 

3. Analysis of avalanches 

A time series X(t) > 0, representing a succession of avalanches with zero quiet times, 
can be constructed numerically from the discrete-time version of Eq. (jlj) by integrating 
the equation using realizations of the fractional Gaussian noise process AW nit)- At 
those times when X(t) drops below zero we consider the avalanche as terminated, and 
a new, independent realization of AWuif) is generated and used to produce the next 
avalanche. From long, stationary time-series generated from the stochastic model and 
from the sandpile model this way, we can construct pdfs V(X) which turn out to give 
almost identical results for the two models (see Fig. [2]). The shape of this pdf is universal 
in the thermodynamic limit: a stretched exponential V(X) ~ exp (-aX'") with fi pa 0.5. 
A different pdf appears if the time-series are constructed by launching the avalanches at 
random times (Poisson-distributed) with characteristic time between launches shorter 

| The BTW model does not exhibit perfect finite-size scaling [5] and hence the scaling xn ~ ^V Dl is 
not valid for very large activity. The effect of imperfect scaling with increasing N can be built into Eq. 0] 
through an TV-dependent drift term. However, the distributions of duration and size of sub-system size 
avalanches (defined by a threshold X c > 0) is not sensitive to this feature of the BWT model. We have 
given a detailed treatment of this problem in [6] . 
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Figure 2. Logarithmic plots of V(X) from simulations of the 2-dimensional BTW 
sandpile for N — 1024. Also shown is V{X) found from simulations of Eq. (j4]), and a 
stretched exponential fit (dashed curve, vertically shifted for visibility). All pdfs are 
scaled to unit variance. 



than the growth time of a system size avalanche. In this case several avalanches may 
run simultaneously, and V(X) from both models are close to the Bramwell-Holdsworth- 
Pinton distribution, which was claimed to be valid for the toppling-activity in the BTW- 
model in [3]. 

Consider a solution of Eq. (j3J) with initial condition X(0) = Y > 0, and let 
P(X, t) be the evolution of the density distribution in X-space of an ensemble of 
realizations of the stochastic process X(t) all launched at activity X = Y at time t — 0. 
Every realization X(t) will sooner or later terminate at a finite time t = r for which 
X(t — 1) > and X(t) < 0, and then we remove it from the ensemble. P(X, t) contains 
information about all commonly considered avalanche characteristics. For example, it 
is easily found from from Eq. (J4j) that, on time scales shorter than the growth time of 
a system-size avalanche, X(t) is a self-similar process with non-stationary increments 
and self-similarity exponent h = 2H [6]. Hence the variance of X(t) with respect to 
P(X,t) will scale as ~ t 2h . That this relation holds for the 2-dimensional BTW model 
can easily be verified through numerical simulation (Fig. 13(a)). 

We can also compute the survival probability p(r) = J °° P(X, r) dX, which is the 
probability that a realization of an avalanche has not terminated at the time r. This 
function is related to the pdf for avalanche durations by Pdur(T) = — p'(r), so that pdur(T) 
is a power law if and only if p(r) is a power law. Fig. [3]^b) shows the function p(r) for 
numerical simulations of the BTW sandpile in the re-scaled coordinates Xn and t^, 
demonstrating that the pdf for avalanche durations does not represent a power law. 
The power-law form p(r) ~ r 5 proposed in [7] can only be obtained as a tangent to the 
log-log plot of p(r) at a given duration time r, and the slope of this tangent depends 
crucially on the duration time r for which this tangent is drawn. 

The situation changes if we let all avalanches terminate when X drops below a 
small threshold X c > as proposed in [1]. In this case avalanche durations are the 
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Figure 3. a) Double- logarithmic plots of the variance of X(t) with respect to the 
pdf P(X,t). The variance grows like t 2h , with h = 2H = 0.74 for times less than 
the duration of a system size avalanche, b) Double-logarithmic plots of the survival 
function p(r) in the re-scaled coordinates Xn and t^, demonstrating that the pdf of 
avalanche durations is not a power-law. The dotted line has slope —0.5. 



return times to the line X = X c , and by changing coordinates to Y = X — X c we 
see that this corresponds to the return times to the time axis of the process given by 
the stochastic differential equation dY(t) = a -J X c + Y(t) dWjjif). For small avalanches 
where X{t)—X c <C X c we can approximate this expression with dY{t) = a y/XZdWH{t), 
i.e. can approximate Y(t) by a fractional Brownian motion with Hurst exponent H. 
Using the result of Ding and Yang [9] on the return times of a fractional Brownian 
motion we get Pdur( r ) ~ t 2 ~ h = r~ L63 . 

Numerical simulations of the BTW model verify this result: The survival function 
p(r) becomes a power law on time scales shorter than a system-size avalanche (see 
Fig. Ufa)), and the slope of the graph in a log-log plot is approximately —0.63, which 
corresponds to a scaling of the pdf for duration times on the form pdur(T) ~ r~ L63 . The 
result is also reproduced by simulations of Eq. (TJJ with an exponentially decaying drift 
term. Fig. H(b) shows the log-log plot of the pdf for duration times in the stochastic 
differential equation and a line with slope —1.63, demonstrating that the avalanche 
statistics in the BTW sandpiles is captured by the stochastic differential equation. 
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Figure 4. a) The survival function for the BTW sandpile in re-scaled coordinates 
Xm and for N = 64, 256, 1024, 2048 when the durations are defined by putting a 
small threshold X c on the toppling activity. The function shows power-law behavior 
with exponent —0.63 for avalanches smaller than system size, b) The pdf for duration 
times from simulations of Eq. (j4|) when avalanches are defined in the same way as for 
the sandpiles. The dotted line has slope —1.63. 



From the scaling p(r) ~ r~ a we can deduce an exponent for the pdf of avalanche size 
as well. On the time scales where the toppling activity can be approximated by a 
fractional Brownian motion Wn{t), the signal disperses with time as X ~ t , the size 



of an avalanche of duration r scales like S(t) ~ J Q t H dt 



~ T 



H+l 



Assuming that the 



pdf for avalanche size is on the form p S i ze (S) ~ S v , the relation p S i ze (5') dS = 
yields T -^( H + 1 )+ H ~ t' 01 ' 1 , so 

H+a+l 2 



(t) dr 



(5) 



H+l H+l 
With H = 0.37 we obtain v = 1.46. 

We also remark that if we omit the drift term and let H = 1/2 and X c = we 
obtain the so-called mean-field theory of sandpiles. In this case the stochastic differential 
equation has a corresponding Fokker-Planck equation 

d -^ = a -^(xp) 

dt 2 dX 2V h 

If we solve this equation on the interval [0, oo) with absorbing boundary conditions in 
X = we can obtain an analytical expression for P(X, t), and from some straightforward 
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algebra we find for large r that Pdur( r ) ~ r ~ 2 [6J- Since X c = we can not approximate 
the toppling activity by a Brownian motion on any scale and thus X(t) diverges like 
~ t h , where h = 2H. By replacing H with h = 2H in Eq. (jSJ) we get p S i ze {S) ~ >S~ 3//2 , 
in agreement with previous mean field approaches [TUJ [TT] . 

4. Concluding remarks 

We point out that the validity of Eq. H]is not restricted to the BTW model. For instance, 
the equation has been verified for the Zhang model [6j[12], though with a different Hurst 
exponent H . Time series of global quantities derived from numerical simulation of 
different sandpile and turbulent fluid systems can be shown to be adequately described 
by Eq. 4, where H and the specific form of the drift term depend on the system at 
hand [6]. 
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